Tutorial 9: Experiments with the BEFWM2

Author

Danet and Becks, based on originals by Delmas and Griffiths

Published

February 28, 2023

The previous tutorial focused on experiments where we manipulated the number of networks and various network parameters. This is one set of things we can change/vary in an in silico experiment. The other set of things we can change are features of the model, such as the shape of the functional response (see Tutorial 7), features of the environment such as the carrying capacity, or even empirical relationships that drive trophic structure and interaction strengths, such as the predator-prey mass ratio.

In this tutorial, we are going to implement three experiments. The first two will be ‘simple’ in that they vary only two things. The final example will implement a large experiment changing five features of the model.

You may want to start a new script in the project. We’ll need the following packages (they are already installed… so we just need using).

using EcologicalNetworksDynamics
using Random, Plots, Distributions, DataFrames, StatsPlots
using EcologicalNetworksPlots, EcologicalNetworks

Experiment 1: Carrying Capacity and the Predator Prey Mass Ratio

Now we are set for our first experiment. Lets first establish the parameters we need to make the food web and do the experiment. We fix S at 20 and C at 0.15. We then create vectors of Z and K.

Z is the predator - prey mass ratio, and defines how much bigger or smaller the predators are from their prey. The data suggest it is between predators are between 10 and 100 times bigger than their prey see Brose et al 2006. This value interacts with setting trophic levels in the model.

The default setting for the models is 1 - i.e. all species are within the same order of magnitude, predators are not bigger than their prey. Here, we creat a vector of values to explore, from predators being smaller, to them being 10 or 100 x larger as the data suggests.

MORE ON HOW THE PPMR interacts with Trophic level?

#Fixed Parameters
S = 20
C = 0.15

# Variable Parameters
Z_levels = [0.1, 1, 10, 100]
K_levels = [0.1, 1, 10, 100]

# run this to get same results as in the document
Random.seed!(123)
TaskLocalRNG()

Now, lets set up the collecting data frame.

df_collect = DataFrame(Z = [], K = [], FinalRichness = [], FinalBiomass = [], FinalStability = [])
0×5 DataFrame
RowZKFinalRichnessFinalBiomassFinalStability
AnyAnyAnyAnyAny

Now, set up the loop to use these variables and generate outputs

for z in Z_levels
    for k in K_levels

        println(" ***> This is iteration with Z = $z and K = $k\n")

        fw = FoodWeb(nichemodel, S; C, Z = z)
        Environment(fw, K = k)

        B0 = rand(S)
        params = ModelParameters(fw)

        out = simulate(params, B0)

        # collect info
        fin_rich = foodweb_richness(out, last = 30)
        fin_bio = total_biomass(out, last = 30)
        stab = population_stability(out, last = 30)

        push!(df_collect, [z, k, fin_rich, fin_bio, stab])
    end
end
 ***> This is iteration with Z = 0.1 and K = 0.1

┌ Info: Species [16] went extinct at time t = 10.333220779245396. 
└ 1 out of 20 species are extinct.
┌ Info: Species [7] went extinct at time t = 58.07510892697318. 
└ 2 out of 20 species are extinct.
┌ Info: Species [8] went extinct at time t = 96.59439608251868. 
└ 3 out of 20 species are extinct.
┌ Info: Species [15] went extinct at time t = 99.81800659688153. 
└ 4 out of 20 species are extinct.
┌ Info: Species [10] went extinct at time t = 173.50017299410476. 
└ 5 out of 20 species are extinct.
 ***> This is iteration with Z = 0.1 and K = 1.0

┌ Info: Species [7, 18] went extinct at time t = 7.79864262878258. 
└ 2 out of 20 species are extinct.
┌ Info: Species [17] went extinct at time t = 12.302588972757128. 
└ 3 out of 20 species are extinct.
┌ Info: Species [14, 19] went extinct at time t = 54.65909534164244. 
└ 5 out of 20 species are extinct.
┌ Info: Species [10] went extinct at time t = 60.53358149873246. 
└ 6 out of 20 species are extinct.
 ***> This is iteration with Z = 0.1 and K = 10.0

┌ Info: Species [5] went extinct at time t = 252.5627381565363. 
└ 7 out of 20 species are extinct.
┌ Info: Species [16] went extinct at time t = 255.6762152446527. 
└ 8 out of 20 species are extinct.
┌ Info: Species [18] went extinct at time t = 5.452845413338652. 
└ 1 out of 20 species are extinct.
┌ Info: Species [17] went extinct at time t = 7.6238324278018. 
└ 2 out of 20 species are extinct.
┌ Info: Species [15] went extinct at time t = 67.45898739065261. 
└ 3 out of 20 species are extinct.
┌ Info: Species [19] went extinct at time t = 70.55144729932225. 
└ 4 out of 20 species are extinct.
┌ Info: Species [11] went extinct at time t = 76.65337831298842. 
└ 5 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 79.93906826954334. 
└ 6 out of 20 species are extinct.
┌ Info: Species [9] went extinct at time t = 206.60040098574206. 
└ 7 out of 20 species are extinct.
 ***> This is iteration with Z = 0.1 and K = 100.0

┌ Info: Species [9] went extinct at time t = 10.740211650168495. 
└ 1 out of 20 species are extinct.
┌ Info: Species [20] went extinct at time t = 21.53600773237581. 
└ 2 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 26.526137967578318. 
└ 3 out of 20 species are extinct.
┌ Info: Species [10, 17] went extinct at time t = 29.65218241510892. 
└ 5 out of 20 species are extinct.
┌ Info: Species [11] went extinct at time t = 40.93973683032451. 
└ 6 out of 20 species are extinct.
┌ Info: Species [7, 3] went extinct at time t = 44.403037167488925. 
└ 8 out of 20 species are extinct.
┌ Info: Species [5, 12] went extinct at time t = 47.422100875122986. 
└ 10 out of 20 species are extinct.
┌ Info: Species [15] went extinct at time t = 53.521946406912726. 
└ 11 out of 20 species are extinct.
┌ Info: Species [13] went extinct at time t = 88.23965733559986. 
└ 12 out of 20 species are extinct.
 ***> This is iteration with Z = 1.0 and K = 0.1

┌ Info: Species [7, 11] went extinct at time t = 33.40288048145466. 
└ 2 out of 20 species are extinct.
┌ Info: Species [16] went extinct at time t = 43.58846848478901. 
└ 3 out of 20 species are extinct.
 ***> This is iteration with Z = 1.0 and K = 1.0

┌ Info: Species [13] went extinct at time t = 71.26258343675524. 
└ 4 out of 20 species are extinct.
┌ Info: Species [9] went extinct at time t = 76.2910085126206. 
└ 5 out of 20 species are extinct.
┌ Info: Species [17] went extinct at time t = 81.47041929142708. 
└ 6 out of 20 species are extinct.
┌ Info: Species [20] went extinct at time t = 250.51901692767748. 
└ 7 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 30.025171952720154. 
└ 1 out of 20 species are extinct.
┌ Info: Species [16, 20, 17, 19] went extinct at time t = 41.4267587744762. 
└ 5 out of 20 species are extinct.
┌ Info: Species [4] went extinct at time t = 61.788208267960805. 
└ 6 out of 20 species are extinct.
┌ Info: Species [11, 18] went extinct at time t = 108.48651859473986. 
└ 8 out of 20 species are extinct.
┌ Info: Species [7] went extinct at time t = 146.3080674035236. 
└ 9 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 296.7203892892406. 
└ 10 out of 20 species are extinct.
 ***> This is iteration with Z = 1.0 and K = 10.0

┌ Info: Species [7] went extinct at time t = 31.97167042395919. 
└ 1 out of 20 species are extinct.
┌ Info: Species [20, 19] went extinct at time t = 44.704826238230446. 
└ 3 out of 20 species are extinct.
 ***> This is iteration with Z = 1.0 and K = 100.0

┌ Info: Species [16] went extinct at time t = 351.1728007870926. 
└ 4 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 93.92059343902136. 
└ 1 out of 20 species are extinct.
┌ Info: Species [13] went extinct at time t = 232.64065602391264. 
└ 2 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 257.55666477477075. 
└ 3 out of 20 species are extinct.
 ***> This is iteration with Z = 10.0 and K = 0.1

 ***> This is iteration with Z = 10.0 and K = 1.0

┌ Info: Species [19] went extinct at time t = 176.22073341629988. 
└ 1 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 447.84231449641817. 
└ 2 out of 20 species are extinct.
 ***> This is iteration with Z = 10.0 and K = 10.0

┌ Info: Species [18] went extinct at time t = 257.2149928441204. 
└ 1 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 362.2028229313007. 
└ 2 out of 20 species are extinct.
 ***> This is iteration with Z = 10.0 and K = 100.0

┌ Info: Species [12] went extinct at time t = 392.3491322903875. 
└ 1 out of 20 species are extinct.
┌ Info: Species [3] went extinct at time t = 108.36651973920425. 
└ 1 out of 20 species are extinct.
┌ Info: Species [5, 4] went extinct at time t = 142.19180110075308. 
└ 3 out of 20 species are extinct.
┌ Info: Species [11, 9] went extinct at time t = 152.16900575546077. 
└ 5 out of 20 species are extinct.
┌ Info: Species [6] went extinct at time t = 162.39218592381178. 
└ 6 out of 20 species are extinct.
┌ Info: Species [15, 10] went extinct at time t = 184.60725874208512. 
└ 8 out of 20 species are extinct.
┌ Info: Species [7] went extinct at time t = 196.60065563305216. 
└ 9 out of 20 species are extinct.
 ***> This is iteration with Z = 100.0 and K = 0.1

┌ Info: Species [16, 13] went extinct at time t = 254.46108282792977. 
└ 11 out of 20 species are extinct.
┌ Info: Species [20] went extinct at time t = 324.694736964769. 
└ 12 out of 20 species are extinct.
┌ Info: Species [19] went extinct at time t = 360.7907924164594. 
└ 13 out of 20 species are extinct.
 ***> This is iteration with Z = 100.0 and K = 1.0

 ***> This is iteration with Z = 100.0 and K = 10.0

 ***> This is iteration with Z = 100.0 and K = 100.0

┌ Info: Species [4] went extinct at time t = 303.7164417912157. 
└ 1 out of 20 species are extinct.

Wonderful. Now we are in a position to learn about two new plotting methods. First, let’s look at the data frame we’ve created.

df_collect
16×5 DataFrame
RowZKFinalRichnessFinalBiomassFinalStability
AnyAnyAnyAnyAny
10.10.115.02.69924-0.0076855
20.11.012.02.13677-0.0011011
30.110.013.02.18971-0.000546612
40.1100.08.01.50359-0.00038262
51.00.114.42.22238-0.136409
61.01.012.13332.64512-0.178709
71.010.016.82.67383-0.259915
81.0100.018.73.54345-0.0968901
910.00.119.44.81281-0.232271
1010.01.019.58.08706-0.101531
1110.010.019.26674.24443-0.0289865
1210.0100.012.66674.70276-0.0880609
13100.00.120.010.8534-0.42463
14100.01.020.041.1869-0.421736
15100.010.019.714.6763-0.750017
16100.0100.020.018.7396-0.433297

Visualising the experiment

One option here is to plot one of our Final Objects as the response variable against the valuse of Z and K. In R, we’d use ggplot2. Here we’ll use StatsPlots as we learned about in Tutorial 5. Can you make this work in the regular Plots syntax?

Let’s first look at a single plot of stability

@df df_collect plot(:K, [:FinalStability], group = :Z, 
    ylabel = "Stabilty", 
    xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line],
    legend = false)

p1 = @df df_collect plot(:K, [:FinalStability], group = :Z, 
    legend = :bottomright,
    ylabel = "Stabilty", 
    xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line])

p2 = @df df_collect plot(:K, [:FinalBiomass], group = :Z, 
    legend = :bottomright,
    ylabel = "Stabilty", 
    xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line])
    
p3 = @df df_collect plot(:K, [:FinalRichness], group = :Z, 
    legend = :bottomright,
    ylabel = "Stabilty", 
    xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line])

# create a layout of 3 graphs stacked on top of each other.
plot(p1, p2, p3, layout=(3,1), legend = false)

Interpretation!